home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / legendre_poly.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  14.4 KB  |  555 lines
  1. /* specfunc/legendre_poly.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_bessel.h>
  26. #include <gsl/gsl_sf_exp.h>
  27. #include <gsl/gsl_sf_gamma.h>
  28. #include <gsl/gsl_sf_log.h>
  29. #include <gsl/gsl_sf_pow_int.h>
  30. #include <gsl/gsl_sf_legendre.h>
  31.  
  32. #include "error.h"
  33.  
  34. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  35.  
  36. int
  37. gsl_sf_legendre_P1_e(double x, gsl_sf_result * result)
  38. {
  39.   /* CHECK_POINTER(result) */
  40.  
  41.   {
  42.     result->val = x;
  43.     result->err = 0.0;
  44.     return GSL_SUCCESS;
  45.   }
  46. }
  47.  
  48. int
  49. gsl_sf_legendre_P2_e(double x, gsl_sf_result * result)
  50. {
  51.   /* CHECK_POINTER(result) */
  52.  
  53.   {
  54.     result->val = 0.5*(3.0*x*x - 1.0);
  55.     result->err = GSL_DBL_EPSILON * (fabs(3.0*x*x) + 1.0);
  56.     return GSL_SUCCESS;
  57.   }
  58. }
  59.  
  60. int
  61. gsl_sf_legendre_P3_e(double x, gsl_sf_result * result)
  62. {
  63.   /* CHECK_POINTER(result) */
  64.  
  65.   {
  66.     result->val = 0.5*x*(5.0*x*x - 3.0);
  67.     result->err = GSL_DBL_EPSILON * (fabs(result->val) + 0.5 * fabs(x) * (fabs(5.0*x*x) + 3.0));
  68.     return GSL_SUCCESS;
  69.   }
  70. }
  71.  
  72.  
  73. int
  74. gsl_sf_legendre_Pl_e(const int l, const double x, gsl_sf_result * result)
  75. { 
  76.   /* CHECK_POINTER(result) */
  77.  
  78.   if(l < 0 || x < -1.0 || x > 1.0) {
  79.     DOMAIN_ERROR(result);
  80.   }
  81.   else if(l == 0) {
  82.     result->val = 1.0;
  83.     result->err = 0.0;
  84.     return GSL_SUCCESS;
  85.   }
  86.   else if(l == 1) {
  87.     result->val = x;
  88.     result->err = 0.0;
  89.     return GSL_SUCCESS;
  90.   }
  91.   else if(l == 2) {
  92.     result->val = 0.5 * (3.0*x*x - 1.0);
  93.     result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
  94.     return GSL_SUCCESS;
  95.   }
  96.   else if(x == 1.0) {
  97.     result->val = 1.0;
  98.     result->err = 0.0;
  99.     return GSL_SUCCESS;
  100.   }
  101.   else if(x == -1.0) {
  102.     result->val = ( GSL_IS_ODD(l) ? -1.0 : 1.0 );
  103.     result->err = 0.0;
  104.     return GSL_SUCCESS;
  105.   }
  106.   else if(l < 100000) {
  107.     /* Compute by upward recurrence on l.
  108.      */
  109.     double p_mm   = 1.0;     /* P_0(x) */
  110.     double p_mmp1 = x;        /* P_1(x) */
  111.     double p_ell = p_mmp1;
  112.     int ell;
  113.  
  114.     for(ell=2; ell <= l; ell++){
  115.       p_ell = (x*(2*ell-1)*p_mmp1 - (ell-1)*p_mm) / ell;
  116.       p_mm = p_mmp1;
  117.       p_mmp1 = p_ell;
  118.     }
  119.  
  120.     result->val = p_ell;
  121.     result->err = (0.5 * ell + 1.0) * GSL_DBL_EPSILON * fabs(p_ell);
  122.     return GSL_SUCCESS;
  123.   }
  124.   else {
  125.     /* Asymptotic expansion.
  126.      * [Olver, p. 473]
  127.      */
  128.     double u  = l + 0.5;
  129.     double th = acos(x);
  130.     gsl_sf_result J0;
  131.     gsl_sf_result Jm1;
  132.     int stat_J0  = gsl_sf_bessel_J0_e(u*th, &J0);
  133.     int stat_Jm1 = gsl_sf_bessel_Jn_e(-1, u*th, &Jm1);
  134.     double pre;
  135.     double B00;
  136.     double c1;
  137.  
  138.     /* B00 = 1/8 (1 - th cot(th) / th^2
  139.      * pre = sqrt(th/sin(th))
  140.      */
  141.     if(th < GSL_ROOT4_DBL_EPSILON) {
  142.       B00 = (1.0 + th*th/15.0)/24.0;
  143.       pre = 1.0 + th*th/12.0;
  144.     }
  145.     else {
  146.       double sin_th = sqrt(1.0 - x*x);
  147.       double cot_th = x / sin_th;
  148.       B00 = 1.0/8.0 * (1.0 - th * cot_th) / (th*th);
  149.       pre = sqrt(th/sin_th);
  150.     }
  151.  
  152.     c1 = th/u * B00;
  153.  
  154.     result->val  = pre * (J0.val + c1 * Jm1.val);
  155.     result->err  = pre * (J0.err + fabs(c1) * Jm1.err);
  156.     result->err += GSL_SQRT_DBL_EPSILON * fabs(result->val);
  157.  
  158.     return GSL_ERROR_SELECT_2(stat_J0, stat_Jm1);
  159.   }
  160. }
  161.  
  162.  
  163. int
  164. gsl_sf_legendre_Pl_array(const int lmax, const double x, double * result_array)
  165. {
  166.   /* CHECK_POINTER(result_array) */
  167.  
  168.   if(lmax < 0 || x < -1.0 || x > 1.0) {
  169.     GSL_ERROR ("domain error", GSL_EDOM);
  170.   }
  171.   else if(lmax == 0) {
  172.     result_array[0] = 1.0;
  173.     return GSL_SUCCESS;
  174.   }
  175.   else if(lmax == 1) {
  176.     result_array[0] = 1.0;
  177.     result_array[1] = x;
  178.     return GSL_SUCCESS;
  179.   }
  180.   else {
  181.     double p_mm   = 1.0;    /* P_0(x) */
  182.     double p_mmp1 = x;        /* P_1(x) */
  183.     double p_ell = p_mmp1;
  184.     int ell;
  185.  
  186.     result_array[0] = 1.0;
  187.     result_array[1] = x;
  188.  
  189.     for(ell=2; ell <= lmax; ell++){
  190.       p_ell = (x*(2*ell-1)*p_mmp1 - (ell-1)*p_mm) / ell;
  191.       p_mm = p_mmp1;
  192.       p_mmp1 = p_ell;
  193.       result_array[ell] = p_ell;
  194.     }
  195.  
  196.     return GSL_SUCCESS;
  197.   }
  198. }
  199.  
  200.  
  201. int
  202. gsl_sf_legendre_Plm_e(const int l, const int m, const double x, gsl_sf_result * result)
  203. {
  204.   /* If l is large and m is large, then we have to worry
  205.    * about overflow. Calculate an approximate exponent which
  206.    * measures the normalization of this thing.
  207.    */
  208.   double dif = l-m;
  209.   double sum = l+m;
  210.   double exp_check = 0.5 * log(2.0*l+1.0) 
  211.                    + 0.5 * dif * (log(dif)-1.0)
  212.                    - 0.5 * sum * (log(sum)-1.0);
  213.  
  214.   /* CHECK_POINTER(result) */
  215.  
  216.   if(m < 0 || l < m || x < -1.0 || x > 1.0) {
  217.     DOMAIN_ERROR(result);
  218.   }
  219.   else if(exp_check < GSL_LOG_DBL_MIN + 10.0){
  220.     /* Bail out. */
  221.     OVERFLOW_ERROR(result);
  222.   }
  223.   else {
  224.     /* Account for the error due to the
  225.      * representation of 1-x.
  226.      */
  227.     const double err_amp = 1.0 / (GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
  228.  
  229.     double p_mm;     /* P_m^m(x) */
  230.     double p_mmp1;   /* P_{m+1}^m(x) */
  231.  
  232.     /* Calculate P_m^m from the analytic result:
  233.      *          P_m^m(x) = (-1)^m (2m-1)!! (1-x^2)^(m/2) , m > 0
  234.      */
  235.     p_mm = 1.0;
  236.     if(m > 0){
  237.       double root_factor = sqrt(1.0-x)*sqrt(1.0+x);
  238.       double fact_coeff = 1.0;
  239.       int i;
  240.       for(i=1; i<=m; i++) {
  241.         p_mm *= -fact_coeff * root_factor;
  242.         fact_coeff += 2.0;
  243.       }
  244.     }
  245.  
  246.     /* Calculate P_{m+1}^m. */
  247.     p_mmp1 = x * (2*m + 1) * p_mm;
  248.  
  249.     if(l == m){
  250.       result->val = p_mm;
  251.       result->err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mm);
  252.       return GSL_SUCCESS;
  253.     }
  254.     else if(l == m + 1) {
  255.       result->val = p_mmp1;
  256.       result->err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mmp1);
  257.       return GSL_SUCCESS;
  258.     }
  259.     else{
  260.       double p_ell = 0.0;
  261.       int ell;
  262.     
  263.       /* Compute P_l^m, l > m+1 by upward recurrence on l. */
  264.       for(ell=m+2; ell <= l; ell++){
  265.         p_ell = (x*(2*ell-1)*p_mmp1 - (ell+m-1)*p_mm) / (ell-m);
  266.         p_mm = p_mmp1;
  267.         p_mmp1 = p_ell;
  268.       }
  269.  
  270.       result->val = p_ell;
  271.       result->err = err_amp * (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(p_ell);
  272.  
  273.       return GSL_SUCCESS;
  274.     }
  275.   }
  276. }
  277.  
  278.  
  279. int
  280. gsl_sf_legendre_Plm_array(const int lmax, const int m, const double x, double * result_array)
  281. {
  282.   /* If l is large and m is large, then we have to worry
  283.    * about overflow. Calculate an approximate exponent which
  284.    * measures the normalization of this thing.
  285.    */
  286.   double dif = lmax-m;
  287.   double sum = lmax+m;
  288.   double exp_check = 0.5 * log(2.0*lmax+1.0) 
  289.                      + 0.5 * dif * (log(dif)-1.0)
  290.                      - 0.5 * sum * (log(sum)-1.0);
  291.  
  292.   /* CHECK_POINTER(result_array) */
  293.  
  294.   if(m < 0 || lmax < m || x < -1.0 || x > 1.0) {
  295.     GSL_ERROR ("error", GSL_EDOM);
  296.   }
  297.   else if(m > 0 && (x == 1.0 || x == -1.0)) {
  298.     int ell;
  299.     for(ell=m; ell<=lmax; ell++) result_array[ell-m] = 0.0;
  300.     return GSL_SUCCESS;
  301.   }
  302.   else if(exp_check < GSL_LOG_DBL_MIN + 10.0){
  303.     /* Bail out.
  304.      */
  305.     GSL_ERROR ("error", GSL_EOVRFLW);
  306.   }
  307.   else {
  308.     double p_mm;                 /* P_m^m(x)     */
  309.     double p_mmp1;               /* P_{m+1}^m(x) */
  310.  
  311.     /* Calculate P_m^m from the analytic result:
  312.      *          P_m^m(x) = (-1)^m (2m-1)!! (1-x^2)^(m/2) , m > 0
  313.      */
  314.     p_mm = 1.0;
  315.     if(m > 0){
  316.       double root_factor = sqrt(1.0-x)*sqrt(1.0+x);
  317.       double fact_coeff = 1.0;
  318.       int i;
  319.       for(i=1; i<=m; i++){
  320.         p_mm *= -fact_coeff * root_factor;
  321.         fact_coeff += 2.0;
  322.       }
  323.     }
  324.  
  325.     /* Calculate P_{m+1}^m. */
  326.     p_mmp1 = x * (2*m + 1) * p_mm;
  327.  
  328.     if(lmax == m){
  329.       result_array[0] = p_mm;
  330.       return GSL_SUCCESS;
  331.     }
  332.     else if(lmax == m + 1) {
  333.       result_array[0] = p_mm;
  334.       result_array[1] = p_mmp1;
  335.       return GSL_SUCCESS;
  336.     }
  337.     else{
  338.       double p_ell;
  339.       int ell;
  340.  
  341.       result_array[0] = p_mm;
  342.       result_array[1] = p_mmp1;
  343.  
  344.       /* Compute P_l^m, l >= m+2, by upward recursion on l. */
  345.       for(ell=m+2; ell <= lmax; ell++){
  346.         p_ell = (x*(2*ell-1)*p_mmp1 - (ell+m-1)*p_mm) / (ell-m);
  347.         p_mm = p_mmp1;
  348.         p_mmp1 = p_ell;
  349.         result_array[ell-m] = p_ell;
  350.       }
  351.  
  352.       return GSL_SUCCESS;
  353.     }
  354.   }
  355. }
  356.  
  357.  
  358. int
  359. gsl_sf_legendre_sphPlm_e(const int l, int m, const double x, gsl_sf_result * result)
  360. {
  361.   /* CHECK_POINTER(result) */
  362.  
  363.   if(m < 0 || l < m || x < -1.0 || x > 1.0) {
  364.     DOMAIN_ERROR(result);
  365.   }
  366.   else if(m == 0) {
  367.     gsl_sf_result P;
  368.     int stat_P = gsl_sf_legendre_Pl_e(l, x, &P);
  369.     double pre = sqrt((2.0*l + 1.0)/(4.0*M_PI));
  370.     result->val  = pre * P.val;
  371.     result->err  = pre * P.err;
  372.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  373.     return stat_P;
  374.   }
  375.   else if(x == 1.0 || x == -1.0) {
  376.     /* m > 0 here */
  377.     result->val = 0.0;
  378.     result->err = 0.0;
  379.     return GSL_SUCCESS;
  380.   }
  381.   else {
  382.     /* m > 0 and |x| < 1 here */
  383.  
  384.     /* Starting value for recursion.
  385.      * Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) ) (-1)^m (1-x^2)^(m/2) / pi^(1/4)
  386.      */
  387.     gsl_sf_result lncirc;
  388.     gsl_sf_result lnpoch;
  389.     double lnpre_val;
  390.     double lnpre_err;
  391.     gsl_sf_result ex_pre;
  392.     double sr;
  393.     const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0);
  394.     const double y_mmp1_factor = x * sqrt(2.0*m + 3.0);
  395.     double y_mm, y_mm_err;
  396.     double y_mmp1;
  397.     gsl_sf_log_1plusx_e(-x*x, &lncirc);
  398.     gsl_sf_lnpoch_e(m, 0.5, &lnpoch);  /* Gamma(m+1/2)/Gamma(m) */
  399.     lnpre_val = -0.25*M_LNPI + 0.5 * (lnpoch.val + m*lncirc.val);
  400.     lnpre_err = 0.25*M_LNPI*GSL_DBL_EPSILON + 0.5 * (lnpoch.err + fabs(m)*lncirc.err);
  401.     gsl_sf_exp_err_e(lnpre_val, lnpre_err, &ex_pre);
  402.     sr    = sqrt((2.0+1.0/m)/(4.0*M_PI));
  403.     y_mm   = sgn * sr * ex_pre.val;
  404.     y_mmp1 = y_mmp1_factor * y_mm;
  405.     y_mm_err  = 2.0 * GSL_DBL_EPSILON * fabs(y_mm) + sr * ex_pre.err;
  406.     y_mm_err *= 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-x));
  407.  
  408.     if(l == m){
  409.       result->val  = y_mm;
  410.       result->err  = y_mm_err;
  411.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(y_mm);
  412.       return GSL_SUCCESS;
  413.     }
  414.     else if(l == m + 1) {
  415.       result->val  = y_mmp1;
  416.       result->err  = fabs(y_mmp1_factor) * y_mm_err;
  417.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(y_mmp1);
  418.       return GSL_SUCCESS;
  419.     }
  420.     else{
  421.       double y_ell = 0.0;
  422.       int ell;
  423.  
  424.       /* Compute Y_l^m, l > m+1, upward recursion on l. */
  425.       for(ell=m+2; ell <= l; ell++){
  426.         const double rat1 = (double)(ell-m)/(double)(ell+m);
  427.     const double rat2 = (ell-m-1.0)/(ell+m-1.0);
  428.         const double factor1 = sqrt(rat1*(2*ell+1)*(2*ell-1));
  429.         const double factor2 = sqrt(rat1*rat2*(2*ell+1)/(2*ell-3));
  430.         y_ell = (x*y_mmp1*factor1 - (ell+m-1)*y_mm*factor2) / (ell-m);
  431.         y_mm   = y_mmp1;
  432.         y_mmp1 = y_ell;
  433.       }
  434.  
  435.       result->val  = y_ell;
  436.       result->err  = (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(y_ell);
  437.       result->err += fabs(y_mm_err/y_mm) * fabs(y_ell);
  438.  
  439.       return GSL_SUCCESS;
  440.     }
  441.   }
  442. }
  443.  
  444.  
  445. int
  446. gsl_sf_legendre_sphPlm_array(const int lmax, int m, const double x, double * result_array)
  447. {
  448.   /* CHECK_POINTER(result_array) */
  449.  
  450.   if(m < 0 || lmax < m || x < -1.0 || x > 1.0) {
  451.     GSL_ERROR ("error", GSL_EDOM);
  452.   }
  453.   else if(m > 0 && (x == 1.0 || x == -1.0)) {
  454.     int ell;
  455.     for(ell=m; ell<=lmax; ell++) result_array[ell-m] = 0.0;
  456.     return GSL_SUCCESS;
  457.   }
  458.   else {
  459.     double y_mm;
  460.     double y_mmp1;
  461.  
  462.     if(m == 0) {
  463.       y_mm   = 0.5/M_SQRTPI;          /* Y00 = 1/sqrt(4pi) */
  464.       y_mmp1 = x * M_SQRT3 * y_mm;
  465.     }
  466.     else {
  467.       /* |x| < 1 here */
  468.  
  469.       gsl_sf_result lncirc;
  470.       gsl_sf_result lnpoch;
  471.       double lnpre;
  472.       const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0);
  473.       gsl_sf_log_1plusx_e(-x*x, &lncirc);
  474.       gsl_sf_lnpoch_e(m, 0.5, &lnpoch);  /* Gamma(m+1/2)/Gamma(m) */
  475.       lnpre = -0.25*M_LNPI + 0.5 * (lnpoch.val + m*lncirc.val);
  476.       y_mm   = sqrt((2.0+1.0/m)/(4.0*M_PI)) * sgn * exp(lnpre);
  477.       y_mmp1 = x * sqrt(2.0*m + 3.0) * y_mm;
  478.     }
  479.  
  480.     if(lmax == m){
  481.       result_array[0] = y_mm;
  482.       return GSL_SUCCESS;
  483.     }
  484.     else if(lmax == m + 1) {
  485.       result_array[0] = y_mm;
  486.       result_array[1] = y_mmp1;
  487.       return GSL_SUCCESS;
  488.     }
  489.     else{
  490.       double y_ell;
  491.       int ell;
  492.  
  493.       result_array[0] = y_mm;
  494.       result_array[1] = y_mmp1;
  495.  
  496.       /* Compute Y_l^m, l > m+1, upward recursion on l. */
  497.       for(ell=m+2; ell <= lmax; ell++){
  498.         const double rat1 = (double)(ell-m)/(double)(ell+m);
  499.     const double rat2 = (ell-m-1.0)/(ell+m-1.0);
  500.         const double factor1 = sqrt(rat1*(2*ell+1)*(2*ell-1));
  501.         const double factor2 = sqrt(rat1*rat2*(2*ell+1)/(2*ell-3));
  502.         y_ell = (x*y_mmp1*factor1 - (ell+m-1)*y_mm*factor2) / (ell-m);
  503.         y_mm   = y_mmp1;
  504.         y_mmp1 = y_ell;
  505.     result_array[ell-m] = y_ell;
  506.       }
  507.     }
  508.  
  509.     return GSL_SUCCESS;
  510.   }
  511. }
  512.  
  513. #ifndef HIDE_INLINE_STATIC
  514. int
  515. gsl_sf_legendre_array_size(const int lmax, const int m)
  516. {
  517.   return lmax-m+1;
  518. }
  519. #endif
  520.  
  521. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  522.  
  523. #include "eval.h"
  524.  
  525. double gsl_sf_legendre_P1(const double x)
  526. {
  527.   EVAL_RESULT(gsl_sf_legendre_P1_e(x, &result));
  528. }
  529.  
  530. double gsl_sf_legendre_P2(const double x)
  531. {
  532.   EVAL_RESULT(gsl_sf_legendre_P2_e(x, &result));
  533. }
  534.  
  535. double gsl_sf_legendre_P3(const double x)
  536. {
  537.   EVAL_RESULT(gsl_sf_legendre_P3_e(x, &result));
  538. }
  539.  
  540. double gsl_sf_legendre_Pl(const int l, const double x)
  541. {
  542.   EVAL_RESULT(gsl_sf_legendre_Pl_e(l, x, &result));
  543. }
  544.  
  545. double gsl_sf_legendre_Plm(const int l, const int m, const double x)
  546. {
  547.   EVAL_RESULT(gsl_sf_legendre_Plm_e(l, m, x, &result));
  548. }
  549.  
  550. double gsl_sf_legendre_sphPlm(const int l, const int m, const double x)
  551. {
  552.   EVAL_RESULT(gsl_sf_legendre_sphPlm_e(l, m, x, &result));
  553. }
  554.  
  555.